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Abstract 

The energy spectrum and the nolinear cascade rates of MHD turbulence is not clearly un- 
derstood. We have addressed this problem using direct numerical simulation and analytical 
calculations. Our numerical simulations indicate that Kolmogorov-like phenomenology with 
^,-5/3 ener gy spectrum, rather than Kraichnan's fc~ 3 / 2 , appears to be applicable in MHD tur- 
bulence. Here, we also construct a self-consistent renomalization group procedure in which the 
mean magnetic field gets renormalized, which in turns yields fc~ 5 / 3 energy spectrum. The nu- 
merical simulations also show that the fluid energy is transferred to magnetic energy. This result 
could shed light on the generation magnetic field as in dynamo mechanism. 
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1 Introduction 

The fluid parcels in a turbulent flow have random motion. However, the random flow velocities 
obeys certain properties. Since most of the talks in the conference dealt with chaos which also yields 
random signals, it is important to contrast the difference between chaos and turbulence. Chaos 
can occurs in a nonlinear system with a few (3-10) degrees of freedom, but a turbulent system has 
many (order millions or more) degrees of freedom. Also, a turbulent system may be chaotic; may 
not be chaotic such as in structures (e.g., vortex street); or it may have both chaos and structure 
coexisting with each other. 

Turbulence is ubiquitous and is of very practical importance. Some of the major applications 
are in mixing, aeroplane and high speed vehicle design, atmospheric flows and weather prediction, 
astrophysical objects like stars, jets, interplanetary medium etc. Even with its wide industrial, prac- 
tical, and theoretical importance, the understanding of turbulence is very weak. There are many 
empirical laws from experiments and simulations. There are some phenomenologies, the most 
famous among these is by Kolmogorov. There only a few mathematically rigorous calculations, pri- 
marily Kraichnan's direct interaction approximation, calculations based on renormalization group 
etc. In this paper we will focus on some of the statistical properties of velocity and magnetic fields 
in a turbulent magnetohydro dynamic (MHD) plasma. We have attempted to review in a concise 
manner some of the phenomenological, numerical, analytical work, and observations from the solar 
wind, with an emphasis on our work (with D. A. Roberts, M. L. Goldstein, J. K. Bhattacharjee, 
V. Eswaran). 

Outline of the paper is as follows. Since the existing MHD turbulence phenomenologies are 
motivated by Kolmogorov's fluid turbulence phenomenology, we review Kolmogorov's arguments in 
section 2. In section 3 we review the existing MHD turbulence phenomenologies. Section 4 contains 
the numerical results, which are compared with the predictions of the turbulence phenomenologies. 
In section 5 we briefly report the cascade rates of fluid and magnetic energies. Section 6 contains 
a renormalization group scheme which provides a self-consistent procedure to obtain the effective 
mean magnetic field. Section 7 contains conclusions. 



where u(x,t) is the velocity field, p(x, t) is the pressure field, p is the density, and v is the viscosity 
of the fluid. We assume that the fluid is incompressible, which implies that V " u ( x >*) = 0- The 
above equation can be made non-dimensional by scaling the variables by large length, velocity, and 
time scales (L, U,T). The ratio of the nonlinear terms (II and III term) and the viscous term (IV 
term) turns out to be UL/v, which is called the Reynolds number. When the Reynolds number 
Re is large (greater that 10,000 or so), we say that the flow is turbulent. Please note that unlike 
chaos, there is no critical Reynolds number above which the flow becomes turbulent; experiments 
exhibit transitions at different Re's. 

Analytic solution to the above equation is not known in turbulent regime. The existence and 
uniqueness of the solution itself is not clear at this moment. From the experiments it is known 



2 Fluid Turbulence 



The fluid flows are described by Navier-Stokes equation, which is 




(1) 
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that the velocity field u(x,i) is random. There is a range of scales present in the flow between 
the energy feeding scale (L) and the dissipation scale (77). This intermediate range is called the 
inertial range. Kolmogorov conjectured that the physics in this range does not depend on the 
energy feeding scale or the dissipation scale; the energy distribution is homogeneous and isotropic; 
the interaction between the velocity modes is local in the Fourier space; and the energy cascade rate 
n(/c) is independent of k under steady state. Using these assumptions, the energy spectrum can be 
deduced immediately using dimensional analysis. We find that one dimensional energy spectrum 
E{k) (J E(k)dk = total energy) is given by 

E{k) = K Ko U 2 ^k- 5 / 3 (2) 

where Kk is an universal constant called Kolmogorov's constant, k is the wavenumber. Note that 
the nonlinear energy cascade II(A;) is equal to the dissipation rate and the energy supply rate of 
the fluid. 

Experiments and simulations fl]] are in reasonable agreement with the above phenomenol- 
ogy. There is a small deviation of the exponent from 5/3 which is attributed to intermittency 
phenomenon. The Kolmogorov phenomenology described above is well supported by the calcu- 
lations based on Direct interaction approximations P], renormalization group techniques ||, |3|, 
self-consistent mode coupling etc. 

Regarding dependence on dimensionality, a brief remark is in order here. In two dimensions, 
there are two different power law exponents, 5/3 for small wavenumber and 3 for large wavenumber 
regimes. The origin of the exponent 3 is because of enstrophy conservation. We will not discuss 
this issue here. 

After this brief introduction to fluid turbulence, now we turn to MHD turbulence. 

3 MHD Turbulence 

The incompressible MHD equations are 

(I T B • V + z T • v) z ± = - - p V Ptot + i/+V 2 z± + v_ v 2 z T 

V • z ± = 

where = u ± b. Here u is the velocity fluctuation, b is the magnetic field fluctuation in velocity 
units (scaled by (47T/9) -1 / 2 ) , Bo is the mean magnetic field in velocity units, and v± = (y ± 
A)/2, where A is the resistivity. The field ptot is sum of thermal and magnetic pressure. 

The Alfven waves are the basic modes of incompressible MHD equations. In absence of the 
nonlinear term (z^ • v) z± > z± are the two independent modes travelling antiparallel and parallel to 
the mean magnetic field. However, when the nonlinear term is present, new modes are generated. 
These modes interact with each other, which results in a turbulent behaviour of the fluctuations. 

Magnetohydro dynamic (MHD) turbulence is more complex than the fluid turbulence. There are 
several MHD turbulence phenomenologies. In the following discussion we will state the arguments 
of Dobrowolny et al. || from which we can obtain the existing MHD turbulence phenomenologies. 
The nonlinear interactions modify the amplitudes of the eddies. We assume that the interaction 
between the fluctuations are local in wavenumber space, and that in one interaction, the eddies zt 
interact with the other eddies of similar sizes for time interval . Then from Eq. (^) , the variation 
Szj^ in the amplitudes of these eddies during these interval is given by 

*4 « -i'Uz k - w 
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In N such interactions, because of their stochastic nature, the amplitude variation will be Az^ 



^N{5z^). Therefore, the number of interactions it takes to obtain a variation equal to its 
initial amplitude z k is 

« \ ? (5) 

(^) 

and the corresponding time T 1 * 1 = iVr^ is 

T± W / !^ + ■ (6) 

The time scale of the energy transfer at wavenumber k is assumed to be T ± . Therefore, the fluxes 
of the fluctuations z k can be estimated as 

^^^(^(^V. (7) 

By choosing different interaction time-scales, one can obtain different energy spectra. Kraichnan 
(1965) and Dobrowolny et al. (1980) || argued that the interacting z k and z k modes will get 
separated because of the mean magnetic field in one Alfven time-scale. Therefore, they chose 
Alfven time scale ta = (kB^)^ 1 as the relevant time-scale and found that 

U + ^U- ^ ^-E + {k)E-(k)k 3 =U. (8) 

If we assume that E + (k) « E~(k), then we immediately obtain 

E+(k) E-(k) ~ (B n) 1/2 AT 3 / 2 (9) 

In absence of mean magnetic field, the magnetic field of the largest eddy was taken as Bq. Kraich- 
nan (1965) also argued that the fluid and magnetic energies are equipartitioned. We refer the above 
phenomenology either as one due to Dobrowolny et al. or the generalized Kraichnan's phenomenol- 
ogy 

If the nonlinear time-scale ~ (kz^) is chosen as the interaction time-scales for the eddies 
z k , we obtain 

& *> (4) 2 (4) *> ( 10 ) 

which in turn leads to 

E ± {k) = K ± (n ± ) 4 / 3 (n T )- 2 / 3 A;- 5 / 3 , (11) 

where If are constants, which we will refer to as Kolmogorov's constants for MHD turbulence. 
Because of its similarity with Kolmogorov's fluid turbulence phenomenology, we refer to this phe- 
nomenology as Kolmogorov-like MHD turbulence phenomenology. This phenomenology was first 
given by Marsch ]7| . It is also a limiting case of a more generalized phenomenology constructed by 
Matthaeus and Zhou, Zhou and Matthaeus ||, which is 

± = A*E+(k)E-(k)k? 
B + ^/kE^k) 



where A is a constant. Here the small wavenumbers (y / kE ± (k) ^> Bo) follow 5/3 spectrum, whereas 
the large wavenumbers (y // kE ± (k) -C Bq) follow 3/2 spectrum. 
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After the discussion on these phenomenologies one would naturally ask which of the phe- 
nomenologies is applicable in the solar wind. From the arguments by Matthaeus and Zhou, and 
Zhou and Matthaeus we see that Kraichnan or Dobrowolny et al.'s phenomenology is expected 
to hold when Bq S> \JkE ± (k); on the other hand Kolmogorov-like phenomenology is expected to 
be applicable when Bq <C y // kE ± (k). We would like to test these scaling arguments from the solar 
wind observations and simulations. In the solar wind, which is a good test ground for MHD turbu- 
lence theories, we find that the exponent of the total energy is 1.69 ± 0.08, whereas the exponent 
of the magnetic energy is 1.73 ± 0.08 ||, somewhat closer to 5/3 than 3/2. This is more surprising 
because Bq ^> \JkE ± {k) for inertial range wavenumbers in the solar wind. Also, in the solar wind 
we do not find a break from 5/3 to 3/2 spectrum as predicted by Matthaeus and Zhou, Zhou 
and Matthaeus H Hence, from the comparison with the solar wind observations, it appears 
that there is some inconsistency in the phenomenological arguments given above. To explore these 
issues further, we performed numerical simulations and some analytic studies. As we will described 
below, the numerical simulations also tend to indicate that the Kolmogorov-like phenomenology, 
rather than Kraichnan's or Dobrowolny et al.'s phenomenology, is probably applicable in MHD 
turbulence. 

We have applied renormalization groups to analyse the above equations. We find that under 
certain assumptions, -Bo gets renormalized as we go from large length-scales to smaller length- 
scales. In other words, Bq appearing in the Kraichnan's or Dobrowolny et al.'s argument must be k 
dependent. This leads to fc -5 / 3 energy spectra, which appears to be consistent with the solar wind 
observations and the simulation results. We will describe these ideas in more detail in the following 
sections. 

Before we proceed further, we point out that the normalized cross helicity a c , defined as {E + — 
E~)/(E + + E~), and the Alfven ratio r^, defined as the ratio of fluid energy and magnetic energy, 
play an important role in MHD turbulence. 

4 Numerical Simulation 

To probe numerically the physics at the inertial range, one usually solves the MHD equations using 
the pseudospectral method with introduction of hyperviscosity and hyperrestivity. For large resolu- 
tion turbulence simulations, spectral method is preferred over finite difference methods because the 
derivatives can be calculated accurately using the Fourier transforms. Note that any error in these 
equations could propagate very fast, which could make the simulation unreliable. Most of our runs 
were performed in two-dimensions because of the expense associated with large three-dimensional 
runs. Fortunately, the scaling arguments described in the earlier section does not depend on the 
dimensionality because the absolute equilibrium theories predict a forward cascade of total energy 
in both two and three dimensions. The simulation method is as follows. 

The equations are time advance in Fourier space. The nonlinear terms are calculated in real 
space, then its Fourier transform is calculated using FFT. For details of the simulation, refer to 
Biskamp and Welter |l0| and Verma et al. [||. Earlier high resolution simulations were performed 
by Biskamp and Welter (1989), Pouquet et al. (1988), and Politano et al. (1989) f(|. In all these 
simulations there was no strong evidence for either 3/2 or 5/3 spectral evidence. 

We have taken v = A, (or = 0) in our simulation. The viscous term has been modified by 
an addition of hyperviscous and hyperresistive term. The modified viscous term is 
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Our v = 5.0 x 10 6 , but k eq = 10.0 which yields Reynolds number of approximately 2 x 10 5 at large 
scales and 500 at small scales. The maximum resolution of our simulation in 2D is 512 2 . The param- 
eter dt = 5.0x 10" 4 . The runs we consider have (B , a c ) = (0, 0), (0, 0.25), (0, 0.9), (1, 0), (1, 0.25), (1, 0.9), (5, 0.25), (, c 
The a c values are for T = 0. We have taken the initial Alfven ratio to be 1.0. The simulation were 
initialized with a box spectrum out to k = 4 for Bq = 0, and a k~ l spectrum out to k = 15 for 
B = 5.0. 

We find that after approximately T = 10 the system has reached a fully developed turbulent 
state, i.e., the energy spectrum shows approximate power laws in the inertial range. We obtained 
the spectral indices by fitting a straight line in the inertial range, and found that it was difficult 
to distinguish between the indices 3/2 and 5/3. Therefore, we decided to study the energy cascade 
rates iP 1 leaving a sphere of radius K in the wavenumber space which is given by 

TS^{K) = -^7 [ K E ± (k, t)dk - 2u [ K [l + (k/k eq )] 2 k 2 E ± (k)dk. (14) 
ot Jo Jo 

The cascade rates H^(K) were calculated numerically and compared with the predictions of the phe- 
nomenologies. The Kolmogorov-like phenomenology predicts that (E^ / E^) / (Jl~ /U + ) 2 = K~/K + , 
while the Dobrowolny et al.'s phenomenology predicts that Il + = In our simulation we find 
that for initial a c = and 0.25, (E^ / E^) / (U~ /H + ) 2 1. For large a c this equality does not hold. 
Note that uncertainties is quite large in our simulations. However, we found in all our runs that 
the larger of E~£ and E^ had larger cascade rate, thus negating Dobrowolny et al.'s predictions. In 
our simulations the mean magnetic field (or the magnetic field of the largest eddies) were at least 
~ 1, whereas the amplitudes of the largest fluctuations in the inertial range (k > 10) were 0.2 
to 0.5 (initial total energy = 1). Hence z^ < Bo, and generalized Kraichnan's model should have 
been applicable, but we find otherwise. Hence we conclude that our numerical result are contrary 
to the predictions of the phenomenologies. The Kolmogorov-like phenomenology appears to hold 
for small a c (with K + = K~) in all -Bo = — 5 regime. However, for large a c , there appears to 
be inconsistency. It is possible that the Kolmogorov-like phenomenology may still hold for large 
a c with constants K + and K~ being unequal and dependent on a c . This issue is under current 
investigation. (The details of the above simulation results are described in Verma et al. ||). 



5 Cascade rates of magnetic and fluid energies 

The origin of magnetic field in the earth, sun, and in the universe is an important problem. It is 
generally believed that the fluid energy gets transferred to the magnetic energy due to nonlinear 
interactions, a mechanism know as dynamo. There are many models to get the desired magnetic 
field configuration in the astrophysical objects. Here we have attempted to investigate the cascade 
rates of magnetic and fluid energies. These studies will shed light on the energy transfer mechanisms 
involved in dynamo mechanism. Since this work is in progress, here we are reporting only the 
preliminary results. 

The equations for fluid and magnetic energies are as follows: 

9 - -.2\ 2 



— - 2ukrj vf = -u • VP - u [u • yu] + u- [b • yb] (15) 
^ - 2Xk 2 ^j b 2 = b- [b • yu] - b- [u • V b] (16) 



The cascade rates by evaluating — L T(k)dk where T(fc)'s are the nonlinear terms appearing in 
the right hand side of the above equations. The pressure term does not yield any energy transfer. 
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The second term in the RHS of Eq. (15) yields the energy transfer from u — ► u (H u ), while the 
third term yields energy transfer u — ► b(H u ^b)- The terms in the RHS of Eq. ( |l6l ) primarily yield 
the energy transfer from b — > u (H t b otal )(tlie transfer from b —>■ b has been ignored). The total energy 
transfer rates coming from u and b spheres of radius K are n^ ota ' and n|°* a ' respectively. We have 
obtained these cascade rates numerically. 

We solve the MHD equations numerically one a 256 x 256 grid using the method described in 
the earlier section. We choose v = A = 10 -5 , k eq = 10, and initial a c = 0,r^ = 1.0 . Note that 
there is no forcing in our simulation. We have carried out our simulation till T = 7.5. On DEC 2000 
4/233 the computer time required to reach T = 7.5 was approximately 4 hours. At T = 7.5 we have 
calculated the cascade rates by calculating the triple correlations. These cascade rates are shown 
in Figure 1. We find that H u ^b is positive, while Tib negative. At /c max we have H u ^b = — n^°* a ' 
Hence, there is a net transfer of energy from the fluid energy to the magnetic energy. The energy 
evolution studies indicate that the total fluid energy (u 2 /2) is decaying, while the magnetic energy 
(b 2 /2) is almost constant in time. Hence, the dissipation rate of magnetic energy appears to be 
balanced by the energy gained from the transfer rate from the fluid energy. It is to be seen when 
the magnetic energy get enhanced. In any case, these simulations show that there is a net transfer 
of energy from velocity field to magnetic field. We also notice that the energy cascade rate Tl u is 
much smaller that both U u ^b and Y[ t b otal . Careful analysis of the cascade rates for various parameter 
range is in progress. 

We have attempted to obtain the perturbative solutions of MHD equations. This is described 
below. 



6 Kolmogorov-like powerlaw using renormalization groups 

There have been earlier attempts of applying renormalization groups Jll] ] and other analytic tech- 
nique, e.g. Eddy-Damped-Quasi-Normal-Markovian approximation, to MHD turbulence. In earlier 
renormalization group calculations, corrections to viscosity and diffusivity are evaluated by the 
nonlinear terms in presence of forcing. In our calculation presented here, we attempt to calculate 
corrections to the mean magnetic field due to the nonlinear terms. 

The basic idea of our calculation is to obtain the effective Bq at higher wavenumbers. In the 
phenomenology of Kraichnan and Dobrowolny et al. 0] the external field or the magnetic field of 
the largest eddy is taken to -Bo- Here we construct a self-consistent scheme in which the effective 
Bq is taken to be the magnetic field of the next-largest eddy. For example, for Alfven waves of 
wavenumber k, the effective magnetic field -Bo will be the magnetic field of the eddy of size k/10 or 
so. This argument is based on the physical intuition that for the scattering of the Alfven waves at 
a wavenumber k, the effects of the magnetic field of the next-largest eddy is much more than that 
of the external field. (For similarity, please note that in WKB method , local inhomogeneity of the 
medium determines the amplitude and phase evolution). We will find in the following discussion 
that the k dependent -Bo yields /c -5 / 3 energy spectra. For simplicity we have taken E + (k) = E~(k) 
and ta = 1- 

The MHD equation in the Fourier space is || 

-«f(k,t)T«(B -k)zf(k,t) = -iM ijm (k) J dp*T(p,t)4(k-p,t) (17) 

where 

kik m 

-^ijm(k) = /cj-Pj m (k); -Pj m (k) = 5i m -p^ , (18) 
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Here we have ignored the viscous terms. The above equation will, in principle, yield an anisotropic 
energy spectra (different spectra along and perpendicular to Bo). Here we modify the above 
equation in the following way to preserve isotropy 



d_ 

dt 



zf (k, t) T i (B k) zf (k, t) = -iM ijm (k) J dpzf(p, t)zf{k - p, t) 



(19) 



This equation can be thought of as an effective equation with an isotropic random mean field. 

The calculation of Bo(k) using the renormalization group is done by a procedure adopted by 
McComb, McComb and Watt, Zhou et al., and others ||. The wavenumber range (ko.-K^), is 
divided logarithmically into N divisions. The nth shell is (k n -i..k n ) where k n = s n k Q {s > 1). Brief 
summary of the RG operation used is given below (For details of the procedure, refer to || and 

@). ' 

1. Decompose the modes into the modes to be eliminated (A; < ) and the modes to be retained 
(fc > ). In the first iteration (ko--k±) = A; < and (ki..kpf) = k > . 



2. We rewrite the Eq. (19) for A; < and /c > . The equation for the retained modes are 



dt 



zf >(k, t) =F i (B (k)k) zf>(k, t) = -iM ijm (k) J dp [zf>(p, t)z±>(k - p, t) 
zj> (p, t)4< (k - p, t) + zj< (p, t)z±> (k - p, t) + zj< (p, t)z±< (k - p, t) 

We get a similar equation for ^^(k, t) modes 



+ 



(20) 



3. The terms given in the second bracket in the RHS of Eq. (|20|) is calculated perturbatively 
and the z i < (k, t) modes are averaged out. The effective equation after this operation is 



izt>(k,t) T i ([Bo(fc) +SB±(k)] k) zt>(k,t) 
= -iM ijm (k) J dp \zj> (p, t)4> (k - p, t) 



(21) 



where 



5B±{k) 



(22) 



The term b2(k,p,q) = kp(l + z 2 )(z + xy),where x,y,z are the cosine of angles between 
(p,q), (q,k), and (k,p) |]. Note that 5B$(k) = 5Bg(k) (for E + = E~ and r A = 1). Let us 
denote 

Bt(k) = B (k) + 6B (k) (23) 

4. We keep eliminating the shells by the above procedure. After n + 1 iterations we obtain 

B n+l (k) = B n (k) + 5B n (k) (24) 

where 

SBJk) - 4 jL*"***" (2?) kski) (25> 

5. We substitute the following forms for E(k) and B n {k) in the Eqs. (|24,25) and solve for B*(k') 
numerically j|, 12]. 

E{k) = aU 2 / 3 k- 5 / 3 
B n (k n k') = a^U^kn^B^k') 
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We start with initial Bo(ki) at each shell, and keep iterating till -B* +1 (A/) ~ B*(k'), that 
is, till the solution converges. We find that B*(k') is approximately proportional to A;'" 1 / 3 . 
Hence, the mean magnetic field scales as /c _1//3 , and the energy spectra scales as /c _5//3 in this 
self-consistent scheme. 

Hence, we see that scaling of Bq leads to A: -5 / 3 energy spectra. This is a self-consistent scheme 
that shows the plausibility of Kolmogorov-like energy spectra for MHD turbulence. Note that 
simulations and the solar wind observations appear to favour the Kolmogorov-like phenomenology. 
Therefore, the above analytical analysis is a promising step toward a proper understanding of the 
statistical theory of MHD turbulence. 



7 Conclusions 

In this paper we have reviewed some of the phenomenological, observational, numerical, and analytic 
work in the statistical theory of MHD turbulence, with emphasis on our studies. We find that the 
solar wind observations and the numerical results are inconsistent with the predictions of the 
existing MHD turbulence phenomenologies. The energy spectrum and the cascade rates appear to 
be closer the predictions of Kolmogorov-like phenomenology even when the mean magnetic field or 
the magnetic field of the largest eddies is large compared to the inertial range velocity and magnetic 
field fluctuations (a region where Kolmogorov-like theory is not expected to hold). 

We have attempted to obtain Kolmogorov-like energy spectrum in MHD turbulence in presence 
of arbitrary Bq by postulating that the effective Bq is scale dependent, unlike what has been taken 
by Kraichnan and Dobrowolny et al. We have constructed a renormalization group scheme and 
shown that the self consistent Bo(k) oc A; -1 / 3 and E{k) oc A;~ 5 / 3 . This analysis has been worked out 
when E + = E~ and = 1. The generalization to arbitrary parameters, or at least to the limiting 
cases, is also planned. We will be able to the get the Kolmogorov's constants for MHD turbulence 
analytically using this procedure; these constants will be useful for the large-eddy-simulations (LES) 
of MHD turbulence. 

We have also carried out a preliminary study of the cascade rates of fluid and magnetic energies. 
We find that there is a net transfer of fluid energy to magnetic energy. We are investigating the 
parameter regimes in which these cascade rates are large enough (more than dissipation rates) so 
that the total magnetic energy increases with time, what is found in dynamo theory. 

Lastly, we would like to mention that the MHD turbulence theories are very important for 
modelling various astrophysical phenomena and plasma processes. We have estimated turbulent 
heating, nonclassical viscosity and resistivity of the solar wind using the MHD turbulence phe- 



nomenologies [ 13 [ . In this light, search for a satisfactory theory of MHD turbulence appears quite 
important. 

We thank all our collaborators, D. A. Roberts, M. L. Goldstein, J. K. Bhattacharjee, and V. 
Eswaran. MKV also thanks V. Subrahmanyam, M. Barma, V. Ravishankar, and D. Sa for numerous 
useful discussions. 
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FIGURE CAPTIONS 



Figure 1 This figure shows the fluxes for a simulation with initial a c = 0. and = 1.0. The fluxes 
are found at t = 7.5 when a c = 0.015 and = 0.4. 
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